      Program Stat
      !*****************************************************************
      include 'accur.for'
      include 'statd.for'
      Parameter(S_mass=1.989d33)
      Parameter(IMax0=1000)
      Dimension r(0:IMax),alpha(0:IMax),theta(0:JMax)
c      real etime,dtime,tarray(2) ! Convex
      Common /phi/phi(0:IMax,0:JMax)
      Common /delta/delta
      Common /pi/pi
      Common /IB/IB/JB/JB/I0/I0
      Dimension phi0(0:IMax0)
      Logical LogIter,LogInit,LogRapid,Logphi

c      i=errtrap('4'x) ! Convex: trap at math. errors
      ! === constants ==>
      LogRapid=.true.
      JB=0 ! spheroidal configurations
      Logphi=.true.
c      Logphi=.false.
      delta=2d0 ! grid parameter
      rho0=2d14
      ! <== constants ===
c      open(1,file='output',access='append')
      open(1,file='output')
c      Call timer(ITime0) ! PC
          write(1,*) '    rho0   ',
     +               'r(IB)**delta',
     +               '     Re    ',
     +               'CM_u/S_mass',
     +               ' CJ_u/1d50 ',
     +               '  T/abs(W) ',
     +               '   gamma_  '
	  write(1,*) '   rhomax  ',
     +               '    phi0   ',
     +               '    phiA   ',
     +               '    phiB   ',
     +               '    CM     ',
     +               '     V     ',
     +               '     CJ    '
	  write(1,*) '     W     ',
     +               '    VT     '
c      close(1)
      pi=acos(-one)
      dr=one/IA
      do i=0,IA
          r(i)=dr*i
      enddo
      dalpha=one/r(IA-1)-one/r(IA)
      do l=0,IMax-IA+1
          alpha(l)=dalpha*l
      enddo
      do i=IA+1,IMax-1
          l=IMax-i
          r(i)=one/alpha(l)
      enddo
      dtheta=half*pi/JMax
      do j=0,JMax
	  theta(j)=dtheta*j
      enddo
c     === initial phi ==>
      if(Logphi) then
          do j=0,JMax
              do i=0,IMax
                  if(i.le.IA) then
                      phi(i,j)=float(i-IA-1)/float(IA)
                  else
                      l=IMax-i
                      phi(i,j)=phi(IA,j)*degree(alpha(l),delta)
                  endif
              enddo
          enddo
      else
          open(16,file='phi',status='old')
          read(16,*) phi0
          close(16)
          do j=0,JMax
              do i=0,IMax
                  if(i.le.IA) then
                      dr0=one/IMax0
                      im=degree(r(i),delta)/dr0
                      ip=im+1
                      omm=im*dr0
                      omp=ip*dr0
                      xm=(omp-degree(r(i),delta))/dr0
                      xp=(degree(r(i),delta)-omm)/dr0
                      phi(i,j)=xm*phi0(im)+xp*phi0(ip)
                  else
                      l=IMax-i
                      phi(i,j)=phi(IA,j)*degree(alpha(l),delta)
                  endif
              enddo
          enddo
      endif
c     <== initial phi ===
      if(JB.eq.0) then ! spheroidal configurations
          I0=0
      else ! thoroidal configurations
          I0=half*(IB+IA)
      endif

      IB=IA
      do while(IB.gt.3) ! ===== IB cycle ====>

c      rho01=1d15
c      rho02=5d14
c      irho0max=20
c      do irho0=0,irho0max ! ===== rho0 cycle ====> 
c          IB=IA
c          rho0=rho01+float(irho0)/max(one,float(irho0max))*(rho02-rho01)

          LogInit=.true.
	  write(*,*) 'rho0, IB', rho0, IB
          C=( phi(IA,JMax)+Fpsi(one)
     +       -phi(I0,JMax)*EOS(zero,3)/EOS(rho0,3)
     +      )
     +      /(one-EOS(zero,3)/EOS(rho0,3))
          phiu=EOS(rho0,3)/(C-phi(I0,JMax))
	  n_max=20
	  LogIter=.true.
	  n=1
          do while(LogIter.and.(n.le.n_max))
	      C0=C
              Call new_phi(r,alpha,theta,dtheta,
     +                     C,phiu,rho0,LogInit,LogRapid)
	      if(dabs((C0-C)/C).lt.1d-7) LogIter=.false.
              write(*,*) n,C,phi(I0,JMax)
	      n=n+1
          enddo
1         Format(i5,4(1pe10.2),1pe28.20)
          Call Integrals(CM,V,CJ,T,W,Eth,Pip,phi,r,theta,C,phiu,
     +                   rho0,omegaA)
	  VT=abs(2d0*T/W+one+3d0*Pip/W)
          rhomax=zero
          do i=0,IA
              do j=0,JMax
		  rhomax=max(rhomax,
     +                     EOS(phiu
     +                         *( C
     +                           -phi(i,j)
     +                           -Fpsi(degree(r(i),
     +                                        delta)*sin(theta(j)))),
     +                          1
     +                     )/rho0)
              enddo
          enddo
c          open(1,file='output',access='append')
c          write(1,*) 'time',etime(tarray),dtime(tarray)
          ! === PC ==>
c          Call timer(ITime)
          IDTime=ITime-ITime0
          ! <== PC ===
          write(1,*) '*** IB, time', IB, IDTime
          G=6.67d-8
          Re=sqrt(phiu/(G*rho0))
          write(1,2) rho0,r(IB)**delta,Re,CM*Re**3*rho0/S_mass,
     +               CJ*sqrt(G)*Re**5*rho0**1.5d0/1d50,
     +               T/abs(W),Pip/Eth+one
	  write(1,2) rhomax,phi(I0,JMax),phi(IA,JMax),phi(IB,JB),
     +               CM,V,CJ
	  write(1,2) W,VT
2         Format(7(1pe11.3))
c          close(1)
c          open(2,file='f.dat',access='append')
          open(2,file='f.dat')
          write(2,2) rho0,r(IB)**delta,Re,CM*Re**3*rho0/S_mass,
     +               CJ*sqrt(G)*Re**5*rho0**1.5d0/1d50,
     +               T/abs(W),Pip/Eth+one
          close(2)
          write(*,*) 'CJ, VT',CJ*sqrt(G)*Re**5*rho0**1.5d0/1d50, VT
          if(IB.le.0.2d0*IA) then
              IB=IB-1
          else
              if(IB.le.0.4d0*IA) then
                  IB=IB-2
              else
                  if(IB.le.0.6d0*IA) then
                      IB=IB-4
                  else
                      IB=IB-8
                  endif
              endif
          endif

      enddo
      do i0=0,IMax0
          om_=float(i0)/IMax0
          im=degree(om_,one/delta)/dr
          ip=im+1
          omm=im*dr
          omp=ip*dr
          xm=(omp-degree(om_,one/delta))/dr
          xp=(degree(om_,one/delta)-omm)/dr
          phi0(i0)=xm*phi(im,JMax)+xp*phi(ip,JMax)
      enddo
***      open(16,file='phi')
***      write(16,*) phi0
***      close(16)
      stop
      end

      double precision Function EOS(x,n)
      !*****************************************************************
      include 'accur.for'
      Parameter(IMax0=86)
      Dimension rho(0:IMax0),p(0:IMax0),gamma(IMax0+1),CK(IMax0+1),
     +          H(0:IMax0),F(0:IMax0)
      Logical LogInit
      Data LogInit/.true./

      if(LogInit) then
          open(3,file='canuto.dat')
          rho(0)=zero
          p(0)=zero
          IMax=0
          do i=1,IMax0
              read(3,*) rho_c,p_c
              if(rho_c.gt.zero) then
                  IMax=IMax+1
                  rho(IMax)=rho_c
                  p(IMax)=p_c
              endif
          enddo
          do i=2,IMax
              gamma(i)=(log(p(i))-log(p(i-1)))
     +                 /(log(rho(i))-log(rho(i-1)))
              CK(i)=exp(
     +            (log(p(i-1))*log(rho(i))-log(p(i))*log(rho(i-1)))
     +            /(log(rho(i))-log(rho(i-1)))
     +              )
          enddo
          gamma(1)=gamma(2)
          gamma(IMax+1)=gamma(IMax)
          CK(1)=CK(2)
          CK(IMax+1)=CK(IMax)
          H(0)=zero
          H(1)=CK(1)*gamma(1)
     +         *(rho(1)**(gamma(1)-one))
     +         /(gamma(1)-one)
          do i=2,IMax
              H(i)= H(i-1)
     +             +CK(i)*gamma(i)
     +              *( rho(i)**(gamma(i)-one)
     +                -rho(i-1)**(gamma(i)-one))
     +              /(gamma(i)-one)
          enddo
          F(0)=zero
          F(1)=CK(1)
     +         *(rho(1)**(gamma(1)-one))
     +         /(gamma(1)-one)
          do i=2,IMax
               F(i)= F(i-1)
     +              +CK(i)
     +               *( rho(i)**(gamma(i)-one)
     +                 -rho(i-1)**(gamma(i)-one))
     +               /(gamma(i)-one)
          enddo
          LogInit=.false.
      endif
      if(n.eq.1) then ! rho(H)
          if(x.gt.zero) then
              i=IMax+1
              do while(x.lt.H(i-1))
                  i=i-1
              enddo
              if(i.eq.1) then
                  EOS=( x*(gamma(i)-one)/(gamma(i)*CK(i))
     +                )**(one/(gamma(i)-one))
              else
                  EOS=( (x-H(i-1))*(gamma(i)-one)/(gamma(i)*CK(i))
     +                 +rho(i-1)**(gamma(i)-one)
     +                )**(one/(gamma(i)-one))
              endif
          else
              EOS=zero
          endif
      endif
      if(n.eq.2) then ! drho_dH(H)
          if(x.gt.zero) then
              i=IMax+1
              do while(x.lt.H(i-1))
                  i=i-1
              enddo
              if(i.eq.1) then
                  EOS=(x*(gamma(i)-one)/(gamma(i)*CK(i))
     +                )**(one/(gamma(i)-one)-one)
     +                /(CK(i)*gamma(i))
              else
                  EOS=( (x-H(i-1))*(gamma(i)-one)/(gamma(i)*CK(i))
     +                 +rho(i-1)**(gamma(i)-one)
     +                )**(one/(gamma(i)-one)-one)
     +                /(CK(i)*gamma(i))
              endif
          else
              EOS=zero
          endif
      endif
      if(n.eq.3) then ! H(rho)
          if(x.gt.zero) then
              i=IMax+1
              do while(x.lt.rho(i-1))
                  i=i-1
              enddo
              if(i.eq.1) then
                  EOS= H(i-1)
     +                +CK(i)*gamma(i)
     +                 *( x**(gamma(i)-one)
     +                  )
     +                 /(gamma(i)-one)
              else
                  EOS= H(i-1)
     +                +CK(i)*gamma(i)
     +                 *( x**(gamma(i)-one)
     +                   -rho(i-1)**(gamma(i)-one))
     +                 /(gamma(i)-one)
              endif
          else
              EOS=zero
          endif
      endif
      if(n.eq.4) then ! p(H)
          if(x.gt.zero) then
              i=IMax+1
              do while(x.lt.H(i-1))
                  i=i-1
              enddo
              if(i.eq.1) then
                  rho_=( x*(gamma(i)-one)/(gamma(i)*CK(i))
     +                 )**(one/(gamma(i)-one))
              else
                  rho_=( (x-H(i-1))*(gamma(i)-one)/(gamma(i)*CK(i))
     +                  +rho(i-1)**(gamma(i)-one)
     +                 )**(one/(gamma(i)-one))
              endif
              EOS=CK(i)*rho_**gamma(i)
          else
              EOS=zero
          endif
      endif
      if(n.eq.5) then ! E(rho)
          if(x.gt.zero) then
              i=IMax+1
              do while(x.lt.rho(i-1))
                  i=i-1
              enddo
              if(i.eq.1) then
                  EOS= F(i-1)
     +                +CK(i)
     +                 *( x**(gamma(i)-one)
     +                  )
     +                 /(gamma(i)-one)
              else
                  EOS= F(i-1)
     +                +CK(i)
     +                 *( x**(gamma(i)-one)
     +                   -rho(i-1)**(gamma(i)-one))
     +                 /(gamma(i)-one)
              endif
              EOS=EOS*x
          else
              EOS=zero
          endif
      endif
      return
      end
